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SENSITIVITY  ANALYSIS 
FOR 

STATIONARY  PROBABILITIES  OF  A  MARKOV  CHAIN 


by 
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University  of  Wisconsin 


ABSTRACT 

This  paper  considers  the  problem  of  evaluating  the  sensitivity  of  a 
steady-state  cost  oc ( 9)  to  underlying  uncertainty  in  a  parameter  vector  0 
governing  the  probabilistic  dynamics  of  the  system  under  consideration.  We 
show  that  the  gradient  Va(0)  plays  a  fundamental  role  in  the  parametric 
statistical  theory  for  Markov  processes.  We  then  survey  numerical  methods 
available  for  evaluating  Va(0)  and  introduce  a  new  Monte  Carlo  estimator 
for  Va(0),  which  is  applicable  to  Markov  processes  of  substantial 
generality. 
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1 .  INTRODUCTION 


Let  X  ■  {X  :  n  >  0}  be  an  irreducible  positive  recurrent  Markov 
n  — 

chain  governed  by  a  transition  kernel  P( 9) ,  where  0  is  a  parameter 
vector  taking  values  in  i^.  If  x(  9)  is  the  stationary  measure  of 
P(9)  and  f(0,x)  is  the  cost  of  running  the  chain  while  in  state  x,  then 
a(0)  =  /f(9,x)  x(9,dx)  is  the  long-run  average  cost  of  running  X  under 
parameter  choice  9.  In  many  applications  settings,  it  is  of  interest  to 
compute  the  sensitivity  of  a  to  (infinitesimal)  changes  in  the  parameter 
9.  Specifically,  it  is  frequently  useful  to  be  able  to  evaluate  Va( 9) , 
the  gradient  of  a(«)  evaluated  at  9  e  1^.  Since  it  is  generally 
impossible  to  analytically  evaluate  Va(9)  (except  for  simple  models),  this 
paper  will  concentrate  on  numerical  methods  for  determining  Va( 0) . 

This  paper  is  organized  as  follows.  In  Section  2,  we  introduce  an 
important  statistical  application  for  these  methods.  We  show  that  the 
numerical  methods  discussed  here  offer  the  opportunity  to  do  statistical 
point,  variance,  and  interval  estimation  for  highly  complex  functionals  of 
analytically  intractable  Markov  processes.  Section  3  is  devoted  to  the 
formal  derivation  of  an  expression  for  7a( 9)  and  describes,  for  finite 
state  Markov  chains,  a  set  of  linear  equations  which  characterizes  Va(9). 
For  complicated  stochastic  processes,  the  corresponding  linear  systems  are 
too  complex  to  solve  via  standard  numerical  methods,  and  Monte  Carlo 
techniques  therefore  become  relevant.  Thus,  Section  4  provides  a  (new) 
Monte  Carlo  estimator  for  7a(9),  which  is  applicable  to  Markov  chains  of 
substantial  generality.  Finally,  Section  5  offers  a  brief  summary  of  the 
paper . 


2.  STATISTICAL  RELEVANCE  OF  THE  GRADIENT 

Suppose  that  the  transition  kernel  P  governing  the  Markov  chain  X 

is  determined  by  a  finite  family  of  distributions  (F^,  . ..,  F^)  =* 

(F, (y , ) ,  F  (y  )),  where  each  F.(y.)  is  a  probability  distribution 

11  mm  i  i  ^ 

associated  with  a  known  parametric  family  in  which  c  *  i.  If 

0  *  (y.,  •••,  Y  )»  then  P  can  be  viewed  as  a  function  of  9,  namely 
1  m 

P  *  P(0). 

In  statistical  contexts,  the  vector  9  £  1**  (d  ■  d ,  +  •••  +  d  ) 

1  m 

is,  in  general,  unknown.  Most  of  the  literature  on  statistical  inference 
for  Markov  processes  has  concentrated  on  estimation  of  the  "true”  parameter 
9*  (i.e.,  estimation  of  9  when  the  observed  chain  X  is  governed  by 
P(9*))  and  on  related  issues  such  as  production  of  variance  estimates  and 
confidence  intervals.  However,  in  many  applications  settings,  it  is  of 
more  practical  importance  to  estimate  not  9*  but  some  associated  steady- 
state  cost  a(0*) . 


(2.1)  EXAMPLE.  Let  X  *  {X  :  n  >  0}  be  the  Markov  chain  consisting  of 

n  — 

waiting  times  of  consecutive  customers  in  the  M/M/l/®  queue.  (See  HEYMAN 
and  SOBEL  (1982)  for  a  description.)  Arrivals  follow  an  exp(Yj)  distri¬ 
bution,  whereas  service  times  are  distributed  exp(Y2)*  Suppose  that  the 
long-run  customer  waiting  time  a( 9)  is  of  importance,  when  9  - 
The  objective  is  to  produce  estimates  for  a(8*),  as  well  as  variance  and 
interval  estimates,  from  observed  inter-arrival  times  Y^,  •••»  as 

well  as  observed  service  times  Y„, ,  ...,  Y„  .  Note  that  in  certain  set- 

21  2  n2 

tings,  the  inter-arrival  times  and  service  times  may  have  been  collected 
from  two  independent  sources,  so  that  no  waiting  times  for  the  system  are 
available.  For  example,  the  queue  might  correspond  to  a  telephone  switch¬ 
ing  system  being  designed,  in  which  historical  inter-arrival  data  exists 
and  service  time  data  for  the  proposed  switching  device  is  available. 
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(2.2)  EXAMPLE.  Virtually  any  general  discrete-event  stochastic  system  can 

be  formulated  as  a  generalized  semi-Markov  process  (GSMP).  A  GSMP  can, 

in  turn,  be  viewed  as  a  Markov  chain  X  ■  (X  :  n  >  0),  where  X  ■  (S  ,C  ) 

n  —  n  n  n 

records  the  "physical  state"  S  (e.g.,  configuration  of  customers  in  a 

n 

queue)  and  clock  readings  C  (e.g.,  remaining  service  times  for  each  of 

^  t  H 

the  customers  in  the  system)  at  the  n  transition  of  the  GSMP.  (For 
further  details,  see  GLYNN  (1983).)  GSMP's  are  characterized  probabilis¬ 
tically  by  certain  distributions  F^,  ...,  F^  governing  the  way  clocks  are 
reset  (e.g.,  service  times  in  a  queue)  and  by  routing  probabilities 
p^,  ...,  p^  (e.g.,  the  proportion  of  customers  who  visit  station  j  after 
receiving  service  at  station  i). 

In  many  applications  environments,  the  distributions  F^,  F^,  ...,  F^ 

and  routing  probabilities  p^,  ...,  p^  are  unknown  and  must  be  estimated 

via  statistical  methods.  If  one  models  the  distributions  F,  ,  ...,  F.  as 

1  l 

belonging  to  parametric  families  (i.e.,  F  ■  F^y^)),  then  the  transition 
function  P  governing  X  can  be  viewed  as  P  ■  P(  0)  ,  where  0  ■ 

(y  ,  ...»  y^ ,  p^,  ...,  p^) .  The  performance  of  a  stochastic  system  is 
often  assessed  by  considering  a  long-run  average  cost  a  for  the  system 
which,  in  this  context,  can  be  regarded  as  a  function  a  -  a( 0)  of  the 
unknown  parameter  0  associated  with  P.  Consequently,  an  important 
statistical  objective  involves  point  and  interval  estimation  of  a(0*), 
where  0*  is  the  "true"  parameter  governing  the  system. 

We  will  now  outline  a  method  for  obtaining  point  and  interval 
estimates  for  a(0*),  which  is  applicable  to  very  general  stochastic 

*  A  A 

systems.  Let  0  *  (0,(n  ),  ...,  0,(n  ))  be  an  estimator  for  0*  - 

11  d  d 

(0*,  ...»  0*)  (n  is  the  sample  size  associated  with  estimation 
1  a  1  ' 

of  0*.)  Such  estimators  0  are  frequently  available  for  complex  systems. 
In  particular,  one  can  often  appeal  to  maximum  likelihood  estimation  (MLE) 
methods  for  estimating  0*.  Under  very  general  conditions,  0  will  be 
asymptotically  normal,  in  the  sense  that 

(2.3)  0  f  N( 0* ,  C(n  ,  ...,  n  ) ) 

1  a 


where  N(9*,  C(n  ,  n  ))  is  a  multivariate  normal  r.v.  with  mean  9* 

1  d  Q  , 

and  covariance  matrix  C(n.  ,  ...»  n  ).  (  ~  denotes  "has  approximately 

1  a 

the  distribution  of".)  In  certain  design  settings  (see  Example  2.1),  the 
data  for  each  of  the  different  components  9*  is  gathered  from  indepen¬ 
dent  sources.  In  this  case,  C(n  ,  ...,  n  )  takes  the  diagonal  form 

1  d 


(2.4) 


C( n  j ,  ...,  n  ^) 


2. 

Vn  i 


2, 

V-d 


If  a  is  continuously  differentiable  in  a  neighborhood  of  9*,  then  a 
Taylor  expansion  of  a  around  9*  shows  that  (2.3)  yields 

(2.5)  o(§)  *  N(a(9*),  7<r(0*)C  C(n. . n.)  7a(0*)) 

i  a 


where  7a(9*)  is  the  (column)  gradient  of  a  evaluated  at  9*.  (This  is 

the  so-called  "delta  method"  of  statistics.) 

Relation  (2.5)  shows  that  if  n, ,  ...,  n.  are  large,  then  ot(  0)  is 

1  .  d 

a  good  point  estimator  for  <*(9*).  Let  C(n^,  ...,  n^)  be  an  estimator 

for  C(n  ,  ...,  n  )  (such  variance  estimators  are  commonly  available  for 
1  a 

MLE  point  estimates  §).  Then,  (2.5)  proves  that 


(2.6) 


v  =  7a(0)t  C(n, ,  ...»  n  )  7a(9) 
l  a 


is  an  estimator  for  the  variance  of  a(9)  and. 


(2.7) 


a  *1/2  *  *1/2 

[a(9)  -  z(  6)v*'  L ,  a<6)  ♦  s(6)v‘“) 


is  an  approximate  10G(1-6)X  confidence  interval  for  a(9*),  where  z(  6) 
is  the  solution  of  p(n(0,1)  z(6)}  ■  1  -  6/2.  Thus,  provided  that  ot(0) 

and  7a( 9)  can  be  evaluated  (either  analytically  or  numerically),  (2.6) 
and  (2.7)  provide  a  solution  to  the  variance  and  interval  estimation 
problems  discussed  above. 

In  the  case  that  the  covariance  matrix  C(n  ,  ...,  n  )  takes  the 

I  d 

A 

form  (2.4),  v  can  be  expressed  as 


(2.8) 


0  ’  J.  faT  a(8))2  • 


Relation  (2.8)  shows  that  the  contribution  of  uncertainty  in  9*  to  the 
variance  of  ct(0)  is  given  by  ( ba( 0)/ 60^) ^  a^/n^.  This  can  be  used  to 
determine  which  component  to  additionally  sample  if  the  current  estimator 
of  a(0*)  is  too  "noisy." 


(2.1)  EXAMPLE  (continued).  Because  of  the  simplicity  of  the  M/M/l/“ 
queue,  a  can  be  analytically  determined  in  closed  form,  namely  a(y  ,y  ) 

■“  J  A  A  “  b 

Y2(Y2-Yi)  for  <  Y2  (•  for  yl  >  Y2).  If  Yj  <  Y2>  (2.8)  reduces 


2  *2/  1 

-  ,«  °i/ni  *  tt-t- 

<VV  (y2-Yj) 


hi  S£/" 2  • 


-2  -2 

where  Oj(o2)  is  a  variance  estimate  for  y^y^)  formed  from 

Y11 . Yini(Y2r  Y2n2}  * 

For  more  complicated  systems,  such  as  that  described  in  Example  2.2, 
a(<)  cannot  be  determined  analytically,  and  so  one  must  turn  to  numerical 
algorithms.  These  algorithms  will  be  described  in  the  remaining  sections 
of  this  paper. 


3.  A  FORMULA  FOR  THE  GRADIENT  OF  THE  STEADY-STATE 

Let  P(9)  be  the  transition  function  for  X  under  parameter  9,  so 

that  P(9,x,A)  is  the  corresponding  conditional  probability  that  X  e  A, 

n+1 

given  that  X  =  x.  For  an  initial  distribution  n( 9) ,  let  P  be  the 

n  y 

probability  measure  on  the  path-space  of  X  associated  with  P(9),  namely 


pA{x„ 

9  0 


£  A. 


O’ 


A  } 
n 


/  n(9,dxQ)  /  P(9,xQ,dxi)  •••  /  P(9,xn_i>dxn) 


A 


If  X  is  Harris  recurrent  under  P(9)  (see  REVUZ  (1984)),  then  there 
exists  a  unique  probability  measure  u( 9)  such  that 


f(9,x)  7t(  9,dx)  P  a.s. 

y 

f(9)'s).  The  measure  x( 9)  is  stationary 

(3.2)  it(9,*)  -  /  P(9,x,»)  it(9,dx)  . 

S 

(S  is  the  state  space  of  X.)  In  fact,  n( 9)  is  the  unique  probability 
measure  satisfying  (3.2).  Our  goal  is  to  numerically  compute  a( 9)  and 
Va(9),  where  a(9)  is  the  steady-state  limit 

(3.3)  a( 9)  -  /  f(9,x)  x(9,dx)  . 

S 

Since  (3.2)  only  determines  x(9)  up  to  a  multiplicative  constant,  it 
is  necessary  to  add  an  additional  constraint  stating  that  the  total  mass 
x(9,S)  equals  1.  The  quantity  a( 9)  is  then  the  unique  solution  of  the 
integral  equation  system 


(3.D  £  I  f(e,x  )  -  / 

k=0  S 

as  n  -*■  ®  (for  a  large  class  of 
for  P(9),  in  the  sense  that 
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(3.4) 


it(0,*)  m  /  P(9,x,»)  tc(  9,dx) 
S 

?t(9,S)  -  1 

a(9)  ■  /  f(9,x)  x(  9,dx)  . 
S 


The  system  (3.4)  is  well  known  and  has  been  extensively  studied.  If  S  is 
finite,  then  P(9)  is  a  finite  matrix  and  (3.4)  becomes 

*(9)C  -  n(  9) C  P(  9) 

(3.5)  Tt(9)Ce  =*  1 

a(9)  =  x( 9) C  f ( 0) 

(all  vectors  are  column  vectors;  e  is  the  vector  consisting  of  l's). 

As  we  shall  see,  a  similar  system  describes  the  gradient  Va(  9)  of 
a .  Let  us  formally  suppose  that  the  transition  function  P( 9)  can  be 
expanded  as 

P( 9  +  hei)  -  P(  9)  +  hQ1(9)  +  o(h) 
th  d 

where  e^  is  the  i  unit  vector  in  I  .  Assume  that  xCSi-he^)  is 
formally  differentiable  at  h  **  0,  so  that  there  exists  a  signed  measure 
t^(9)  such  that 

(3.6)  u(9  +  he^)  -  n(  9)  +  hr^  9)  +  o(h)  . 

The  stationarity  equation  (3.2)  implies  that  r^(  9)  must  satisfy 

(3.7)  n1(91,dx)  -  /  ^(9, dx)  P(  9,x,  •)  -  /  Q^e.x.O  u(  9,dx) 
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JVV  ■  ' 


-4 


A.V.V.V.-.lV 


•vV/VV  >V' 


v'i 


(formally  differentiate  both  sides  of  (3.2)).  (The  equation  (3.7)  is 
Poisson's  equation  for  the  kernel  P(9).)  These  formal  calculations  can  be 
made  rigorous,  even  in  general  state  space;  such  arguments  will  appear 
elsewhere. 

In  finite  state  space,  the  arguments  are  more  straightforward  and  have 

previously  appeared  in  SCHWEITZER  (1968),  GOLUB  and  MEYER  (1986),  and  MEYER 

and  STEWART  (1986).  We  give  a  very  elementary  proof  in  the  Appendix  to 

this  paper;  our  argument  uses  only  elementary  Markov  chain  theory.  Note 

that  in  finite  state  space,  (3.7)  becomes  q  ( 9)C(I-P(  9))  =  Q^(9). 

This  does  not  uniquely  identify  p  (9),  since  p  ( 9)  +  6x( 9)  also 

^  1  t 

satisfies  the  equation,  for  all  6.  Note  that  since  it(  9)  e  *  1  for  all 
9,  it  follows  that  r|^(9)te  *  0  (see  (3.6)).  Let  II(  9)  be  the  matrix  in 
which  all  rows  are  identical  to  n(9).  It  is  easily  verified  that  since 
ni(9)te  =  0,  ri^( 9) C  11(9)  *  0.  Consequently,  t)^(  9)  also  satisfies 

(3.8)  n  (0)t(I  -  PO)  +  n(9))  =  n(9)t  9)  . 

It  is  well  known  (see  KEMENY  and  SNELL  (1960),  p.  100)  that  ( I-P(  9)+II(  9) ) 
has  an  inverse,  called  the  fundamental  matrix,  which  we  shall  denote  F( 9) . 
Hence,  in  finite  state  space,  the  ith  component  of  Va( 9)  can  be 
computed  as  the  solution  of  the  system 

h,(e)t  -  •rc(9)t  Q  (0)  F(  9) 

(3.9)  *  t  t 

^-a(9)  -  x(9)C  f'(9)  +  Ti1(0)t  f (  9) 

where  f^(9)  is  the  vector  in  which  the  jth  component  is  ^(Q.j)/^©^* 
Consequently,  when  S  is  finite,  the  systems  of  linear  equations 
(3.5)  and  (3.9)  may  be  solved  numerically  to  obtain  a( 9)  and  Va( 9) .  If 
S  is  not  finite  (or  if  the  number  of  elements  in  S  is  large),  numerical 
methods  not  dependent  on  explicit  solution  of  linear  equations  must  be 
considered.  In  the  next  section,  we  show  how  Monte  Carlo  methods  can  be 
used  to  advantage  here. 
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4. 


MONTE  CARLO  EVALUATION  OF  STEADY-STATE  GRADIENTS 

A  critical  assumption  underlying  the  analysis  of  this  section  is  that 
it  is  possible  to  generate  sample  trajectories  of  X  under  the  measure 
Pg .  For  the  examples  that  we  have  in  mind  (see  particularly  Example  2.2), 
this  assumption  is  clearly  in  force. 

Assuming  now  that  X  has  distribution  P  ,  relation  (3.1)  states 

that 

n-1 

(4.1)  -  l  f(9,X,  )  -  a(9)  P.  a.s. 

n  k=0  k  9 


as  n  •*  “.  In  other  words,  rather  than  solving  the  integral  equation 
system  (3.4),  one  may  numerically  approximate  a( 0)  by  the  sample  average 
appearing  on  the  left-hand  side  of  (4.1).  The  simplicity  of  this  numerical 
procedure,  as  well  as  its  broad  applicability,  is  the  source  of  the  power 
of  the  Monte  Carlo  method.  Our  objective  here  is  to  obtain  a  similar  Monte 
Carlo  algorithm  for  evaluation  of  the  gradient  Va( 0) . 

Observe  that  (at  least  formally)  we  have 

(4.2)  -£§-  a( 0)  -  /  £§-  f(0,x)  *(  0,dx)  +  /  f(  0,x)  1^(0, dx)  . 

i  S  i  S 

A  Monte  Carlo  estimator  for  the  first  term  appearing  on  the  right-hand  side 

(4.2)  is  given  by  the  sample  mean 


(4.3) 


1  n_1  a, 

z  l  -si-  • 


“  k:0  sei 


It  remains  to  obtain  an  estimator  for  the  second  term. 


As  in  the  finite  state  space  context,  one  expects  that  the  signed 
measure  r) ^( 0 )  will  satisfy  ti^(0,S)  -0.  As  a  consequence,  it  follows 
from  (3.7)  that  p  (0)  should  satisfy 


(4.4) 


T|  C © , * )  -  /  h  (0,dx)  P(  9,x ,  •)  +  /  r\  (9,dx)  it(  9,  •) 

1  s  1  s 

=  /  Q  (8,x,*)  n(  e,dx)  . 

S  1 

Letting  n(9)  be  the  operator  11(8, x,*)  *  7t(0,*),  one  can  write  (4.4) 
symbolically  as 

(4.5)  n  (0)(I  -  P(9)+  n(  0))  =  n(6)  Q.(0)  . 

(This  is  the  general  state  space  analogue  of  (3.8).)  The  formal  inverse  of 
(I-P(0)+n(0))  is  given  by 


l  (p(0)  -  n(8))k  . 

k=0 

Because  of  the  stationarity  of  tc(  0)  and  the  independence  of  Il(G,x,*) 

k  k 

from  x,  it  follows  that  ( P(  0)-Il(  0) )  =  P(  0)  -  n(  9)  ,  for  k  >_  1 .  Hence, 

a  formal  analysis  of  (4.5)  shows  that 


ri  ( 0)  =  n(0)  Q  ( 9)  +  l  u(0)  Q  (0)(pk(0)  -  n(0))  . 

1  1  k=l  1 

For  the  same  reason  that  ti^(9,S)  m  0,  Q^(9,x,S)  *  0  and  hence 
Q  (6)il(0)  =  0.  Consequently, 

(4.6)  t).  ( 0 )  f  (  0)  *  l  n(0)  Q  (0)  P(0)K  f(  0)  . 

1  k=0  X 

Suppose  that  the  measures  P(»,x,dy)  are  absolutely  continuous  with 
respect  to  P(9,x,dy)  in  a  neighborhood  of  0.  Then,  one  expects  that 
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Q^(9,x,dy)  has  a  density  with  respect  to  P(9,x,dy),  call  it  q^(0,x,y).  A 


typical  term  on  the  right-hand  side  of  (4.6)  then  takes  the  form 


/  x(0,dx)  /  q,(9,x,y)  P(9,x,dy)  /  P  (0,y,dz)  f(0,z) 
S  S  S 


which  can  be  represented  probabilistically  as  an  expectation: 


E0[qi(0.Xo,Xl)f(0,Xk+1)] 


where  Efl(*)  is  the  expectation  corresponding  to  PQ,  and  P.  is  the 


probability  on  path-space  associated  with  initial  distribution  x( 9)  and 
transition  function  P(9).  Thus,  the  second  term  in  (4.2)  has  the  formal 


representation 


(4.7) 


I  e„tq«>,x  x  )«e,x  >1  . 

k-0  °  1  *  1 


The  formula  (4.7)  is  the  key  to  the  Monte  Carlo  analysis. 

Each  term  in  (4.7)  can  be  consistently  estimated  (under  suitable 
hypotheses)  via 


when  X  evolves  according  to  transition  function  P( 0)  (regardless  of  X's 
initial  distribution).  In  order  to  estimate  the  infinite  sum,  a  standard 


device  is  to  consider  an  estimator  of  the  form 


(4.8) 


£(n)  n-Jl(n)-l 


k-o  n'A(n)  j-o  1 J  J+1 


j+k+l ‘ 


s,  mmvm 


where  the  truncation  point  l( n)  is  keyed  to  the  sample  size  n  in  such  a 
way  that  SL(n)  ®  with  Jl(n)/n  >  0.  The  particular  choice  of  Jl(n) 
effects  a  compromise  between  bias  and  variance  effects  in  estimating  the 
infinite  sum  (4.7). 

Since  q^ ( q ,x , y)  is  generally  easily  computable  (for  S  countable, 
qi(8.j.k)  -  (5P(0,j  ,k)/ae^)  •  p(9,j,k)  S,  (4.7)  provides  a  Monte  Carlo 
solution  to  estimating  the  appropriate  gradient. 

It  turns  out  that  (4.6)  is  closely  related  to  a  formula  which  one 
obtains  when  one  uses  likelihood  ratio  change-of-measure  ideas  to  evaluate 
gradients.  These  connections  will  be  explored  more  fully  in  a  future 
paper . 


5 .  SUMMARY 

We  have  shown  that  the  gradient  Va( 9)  of  steady-state  quantity  a 
plays  a  critical  role  in  the  variance  and  interval  estimation  theory  for 
steady-state  estimators  ct(9)  of  complex  stochastic  systems.  In  some 
sense,  the  large-sample  variance  and  interval  estimation  theory  is  fully 

A  A 

solved  given  that  one  can  evaluate  a(9)  and  Va(0).  Numerical  methods 
for  dealing  with  a(9)  when  the  system  is  Markov  are,  of  course,  well 

A 

known.  However,  numerical  algorithms  for  evaluating  Voc(9)  are  a  recent 
development.  We  have  therefore  provided  a  self-contained  exposition  of  the 
relevant  theory,  and  discuss  both  Monte  Carlo  (see  (4.3)  and  (4.8))  and 
non-Monte  Carlo  (see  (3.9))  approaches  to  solving  the  problem. 

APPENDIX 

Let  P(*)  be  a  family  of  n  *  n  stochastic  matrices  which  are: 

(i)  irreducible  in  a  neighborhood  of  9 
(ii)  differentiable  at  9. 

Under  (I),  P( • )  has  a  unique  stationary  distribution  n( •)  in  a 
neighborhood  of  9.  Our  goal  is  to  rigorously  verify  the  first  equation  in 

(3.9) . 

Given  the  existence  of  the  inverse  matrix  F(9)  *  (I-P(  9)+Il(  9)) 

(3.9)  follows  immediately  once  the  differentiability  of  x( 9)  is 
established.  Note  that  for  h  sufficiently  small, 

n(9+hei)  -  u(0)  =  x(  9+hei)t[p(  0)  +  hQi(9)  +  o(h)]  -  x(  9)t  p(  9) 
so 

(x(9+he^)  -  n(0)]t(i-p(9))  -  hu(9+he1)t  Q  (0)  +  o(h) 

(note  that  o(h)n(9+he^)  -  o(h)  since  all  terms  in  i^&fhe^)  are  uniform¬ 
ly  (in  h)  bounded  by  1).  Since  FI(9)  has  identical  rows  and  itfD+he^) 
is  stochastic  for  h  ^  0,  it  follows  that  [  n(  9fhe^)-x(  9)]  IT(  9)  ■  0.  Hence, 


(Al)  [n(0+he^)  -  n(9)]t  -  hi^&fhe^*  Q±(  ©)  F(  e)  +  o(h)  . 

Again,  since  n(9+he^)  is  uniformly  bounded  in  h,  it  is  evident  from  (Al) 
that  u(0+he^)  is  continuous  at  h  *>  0.  Thus,  (Al)  implies  that 

[n(9+he  )  -  u^)]*  -  hit(0)C  0)  F(  9)  +  o(h) 
i.e.,  rli(0)t  -  n(0)C  Qi(0)  F(  0)  , 

which  is  the  required  result. 
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